Finite-element analysis of contact between elastic self-affine surfaces 
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Finite element methods are used to study non-adhesive, frictionless contact between elastic solids 
with self-affine surfaces. We find that the total contact area rises linearly with load at small loads. 
The mean pressure in the contact regions is independent of load and proportional to the rms slope 
of the surface. The constant of proportionality is nearly independent of Poisson ratio and roughness 
exponent and lies between previous analytic predictions. The contact morphology is also analyzed. 
Connected contact regions have a fractal area and perimeter. The probability of finding a cluster 
of area a c drops as a~ T where r increases with decreasing roughness exponent. The distribution of 
pressures shows an exponential tail that is also found in many jammed systems. These results are 
contrasted to simpler models and experiment. 

PACS numbers: 



I. INTRODUCTION 

The forces of friction and adhesion between two sur- 
faces are determined by the interactions between atoms 
at their interface. Most surfaces are rough enough 
that atoms are only close enough to interact strongly 
in areas where peaks or asperities on opposing sur- 
faces overlap. Experiments 

B II H II and theory 

HEHBIllIiinilll show that this real area of con- 
tact A is often much smaller than the projected area Aq 
of the surfaces. They have also correlated 0,0, El the in- 
crease in friction with normal load W to a corresponding 
increase in A. 

Given the importance of A it is not surprising that 
there have been many theoretical studies of the factors 
that determine it in different limits. Some have exam- 
ined the plastic limit where the local pressure is large 
enough to flatten asperities The mean pressure in 
the contacts is W/A and in the simplest model this has a 
constant value that is proportional to the hardness. The 
resulting linear relation between W and A is often given 
as an explanation for the linear rise of friction with load 

urn 

The behavior of elastic contacts is more complicated. 
In the Hertzian limit where friction and adhesion are ig- 
nored, the contact area between a sphere and a flat rises 
only as W 2 ' 3 . However, Greenwood and Williamson's 
pioneering work |l4j showed that a nearly linear rela- 
tion between W and A was obtained if one considered 
a large number of asperities of different height. While 
this calculation assumed spherical asperities with a uni- 
form radius, a linear relation was also obtained when 
it was extended by Bush et al. 15] to include both a 
distribution of radii and aspherical asperities. These cal- 
culations consider explicit probability distributions for 
peaks and sum the Hertzian contact areas calculated for 
each peak without including correlations between peaks. 
Persson Q has recently presented a different approach, 
motivated by the fact that many surfaces have roughness 
on all length scales that can be described by self-affine 



scaling 0, 0, 0] . His calculation considers the scaling 
of stress and contact area with the lengthscale to which 
surface features are resolved. Nevertheless, his final re- 
sult for A only differs from the earlier work of Bush et 
al. ^3 hy a constant factor of it/ 2. 

All of the above theories treat correlations between 
contacting regions approximately. Models based on as- 
perities ignore spatial correlations in asperity heights 
fpH HH . Persson's calculation includes height correla- 
tions and becomes exact in the limit of complete contact 
0. However, when A < Aq, correlations between local 
pressures and contacts are only included in an average 
way. As a result, screening of small bumps on the side 
of larger bumps, may not be included completely. Such 
correlations could change A, and also the distribution of 
pressure along the surface. The range over which A is 
proportional to load is also difficult to determine from 
analytical theories. Indeed, as discussed below, the work 
of Bush et al. 0| seems to suggest that the linear region 
is confined to infinitesimally small loads for a self-affine 
surface. 

A recent numerical study by Borri-Brunetto et al. 
calculated A and the spatial distribution of contact ar- 
eas using a method that is restricted to the limit of zero 
Poisson ratio. They found a substantial range where A 
rose linearly with load, but did not test analytic results 
for the slope specifically. As in Persson's work their 
focus was on the change in results with increasing res- 
olution of surface roughness. Batrouni et aZ.|lOj| have 
considered the scaling of contact area with load using a 
similar method. They find a slight deviation from linear- 
ity, but rule out larger deviations predicted by previous 
scaling arguments 

In this paper, we present numerical calculations of con- 
tact area and pressure distributions for a wide range of 
Poisson ratios, loads, system sizes, self-affine scaling ex- 
ponents and roughness amplitudes. We first show that A 
has a well-defined thermodynamic limit, that is for fixed 
small scale roughness the fraction of area in contact at a 
given average normal pressure is independent of system 
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size. The ratio of A over load is shown to scale as the 
inverse of the rms surface slope, with a coefficient that 
lies between the results of Bush et al. 0] and Persson 
0. The dependence of A/W on roughness exponent and 
Poisson ratio is also obtained. This allows A to be pre- 
dicted for any elastic rough surface. 

We next describe the contact morphology and distri- 
bution of connected contact areas. The probability of 
finding a connected region of area a c falls off as a power 
law, P(a c ) cx a~ T , where r depends only on roughness 
exponent H . For H < 0.9, r is greater than 2, and the 
mean contact size is always comparable to the resolu- 
tion of the calculation. As first noted by Greenwood and 
Williamson the linear rise in A with load reflects 
a linear increase in the number of contacts without any 
increase in their mean size or probability distribution. 
Our results are contrasted to a common model where 
contacts form in the regions where undcformcd surfaces 
would overlap. This approximation has been used to in- 
terpret optical images of contacts 0, Q > an d by Green- 
wood and Wu as a method of estimating the statistics 
of asperity sizes 0- We show that it gives much too 
large a total contact area, and a qualitatively different 
distribution of a c with r < 2. Optical methods may in- 
clude regions that are merely close to touching as part of 
the contact. We find that this can explain the observed 
discrepancy between experiment and calculation. 

The distribution of contact pressures p is also studied. 
We find that results for all system sizes, roughness ampli- 
tudes, and roughness exponents collapse onto a universal 
curve when the pressure is normalized by its mean value. 
The mean value, W/A, can be obtained from the rough- 
ness amplitude as described above. The distribution de- 
creases monotonically with increasing p, and has an expo- 
nential tail at large p. In contrast, approximate analytic 
results for the pressure distribution decay as a Gaussian 
at large pressures Possible connections to exponen- 
tial stress distributions in jammed systems |2(], l2ll 12*2^ 
are mentioned. 

In Sec. II, we describe the numerical procedures used 
to generate self-affine surfaces, mesh them, and deter- 
mine their deformation using an explicit dynamic finite- 
element method. In Section IIIII we present numerical 
results from the contact analysis, including the contact 
area, contact morphology and pressure distribution. The 
final section presents conclusions from the current study 
and discusses avenues for future research. 



II. NUMERICAL SIMULATION 

A. Geometry 

As in previous analytical and numerical calculations 
we use a well-known result from contact mechanics to 
simplify the geometry. If there is no friction or adhe- 
sion between two rough surfaces and the surface slope is 
small, then elastic contact between them can be mapped 



to contact between a single rough surface and a rigid flat 
plane [l3|. The effective modulus E' that controls the 
contact area is given by 1/E' = (1 — vf)/Ei + (l — v%)/E2 
where i>i and Ei are the Poisson ratios and Young's mod- 
uli of the two surfaces. The height h of the new rough sur- 
face is given by the difference between the local heights of 
the original undeformed surfaces. We consider the case 
where both of the surfaces have the same self-affine scal- 
ing properties but are uncorrelated. Then the height h 
has the same scaling properties, but with a larger ampli- 
tude. 

Many surfaces have roughness on all length scales that 
can be described by self-affine fractal scaling. Unlike self- 
similar fractals, self-affine fractal surfaces exhibit differ- 
ent scaling normal to the interface than along it. We 
consider geometries with periodic boundary conditions 
in the x — y plane and specify the surface by the height 
h along the z-axis. For a self-affine surface, variations in 
the height over a lateral length scale £ rise as £ H where 
H < 1 is called the Hurst or roughness exponent. Since 
H < 1, the surface looks smoother at larger length scales. 
Many researchers also specify the scaling by an effective 
fractal dimension d—H, where d is the spatial dimension. 

To generate three-dimensional self-affine fractal sur- 
faces with rms roughness A at small scales, we adopted 
the successive random midpoint algorithm of Voss [iH 
l23l |. The x — y plane is divided into a uniform square 
grid with unit spacing and L nodes along each axis. At 
the first step, the center of the entire grid is displaced by 
a height chosen at random from a Gaussian distribution 
of width £ H A, where £ = is the distance of the 

center from the corners. This center point then becomes 
one corner of new squares that are rotated by 45° and 
have a new corner to center distance £ that is smaller by 
a factor of y/2. The center of each new square is assigned 
a height equal to the average of the corner heights plus a 
random number chosen from a Gaussian of width £ H A. 
This process is iterated down to £ = 1, guaranteeing that 
the variation in height scales with £ in the appropriate 
manner. Figure ^ shows a typical self-affine surface with 
H = 0.5. Note that the height variation is enhanced by 
a factor of ten to make it visible in the figure. 

It is common to test the scaling properties of self-affine 
surfaces by calculating the Fourier transform C(q) 

C(q) = ( 2 *")~ 2 / rf 2 rexp[-iq-r] C(r) (1) 

of the height-height correlation function 

C(r) =< h(r + r')h(r') > (2) 

where (h(r')) — 0, r and r' are vectors in the x — y 
plane, and the brackets indicate an average over r'. The 
function C(q) should be isotropic and decay as q~ 2 ( 1+H \ 
We have generated surfaces with H between 0.3 and 0.9 
and verified that the corresponding C(q) has the correct 
power law scaling. For a given A, the values of C(q) at 
large q are very insensitive to the random seed. How- 
ever, there are large fluctuations at small q where the 



3 




FIG. 1: A self-afnne fractal surface image (256x256) generated 
by the successive random midpoint algorithm. Heights are 
magnified by a factor of 10 to make the roughness visible, and 
the color varies from dark (blue) to light (red) with increasing 
height. 

roughness is dominated by the first few random numbers 
that are chosen at the largest length scales. These first 
random numbers also dominate the mean-squared height 
variation: 

C(r = 0)=<|Mr')| 2 >, (3) 

and we find large variations in C(0) for surfaces with 
the same small scale roughness A. As discussed below, 
the contact area is determined by A and is relatively 
insensitive to fluctuations in C(0). 

B. Mesh Generation 

As noted above, numerical simulations are done for a 
rough elastic surface contacting a perfectly rigid flat sur- 
face. A typical finite element mesh is illustrated in Figure 
13 The mesh is discretized with ten-node tetrahedral el- 
ements. These elements contain three integration points 
and quadratically interpolate the displacement field. A 
coarse mesh is used for the rigid surface to improve nu- 
merical efficiency. A fine mesh for the elastic surface is 
prepared in two stages. 

First, a fine mesh for a flat surface is obtained using a 
longest edge propagation path refinement scheme, which 
is ideally suited to obtain strong mesh gradations as well 
as preserving a high mesh quality [24| . A cube of side L = 
2™ is initially filled with a coarse mesh. Each tetrahedral 
element at the outer surface is then divided to produce 
twice as many surface nodes and the mesh is refined. This 
process is repeated until surface nodes form a uniform 
square grid of unit spacing. Using this technique, meshes 




FIG. 2: Geometry of a finite element mesh in an elastic body 
(top) with a self-affine surface that is pushed down on a flat, 
rigid substrate. Periodic boundary conditions are applied in 
the plane of the interface. 

with L up to 512 are created. The resulting grid contains 
512 x 512 surface nodes, about 911,000 total nodes and 
568, 000 elements. 

Next, the desired surface heights h(x,y) are imposed 
onto the contact surface. Moving only the surface nodes 
produces badly distorted elements, that would at best 
require impractically small time steps and at worst pro- 
duce negative Jacobians. Thus all nodes are moved by 
a fraction of the local height that depends on the initial 
height zq of the node above the bottom of the elastic cube 
(Fig. |2J. The magnitude of the change, Az, decreases 
to zero at the top of the cube so that the top surface 
remains flat. The specific form for the displacement is 
Az(x, y, Zq) = h(x,y) * (L — zq) 11 where a = 6 usually 
gives good meshes. 

C. Finite Element Simulation 

The goal is to determine the equilibrium contact ge- 
ometry at a given load. An implicit approach is too 
memory intensive for the system sizes of interest. In- 
stead we use an explicit integration algorithm combined 
with a dynamic relaxation scheme. Three different algo- 
rithms were compared to insure accuracy. In the first, 
the top surface is given a small velocity and its impact 
with the bottom surface is followed. In the second, the 
displacement of the nodes at the top of the elastic cube 
(Fig. [2J is incremented at a fixed rate or in small discrete 
steps. In the third, a constant force is applied to each 
of these nodes and gradually incremented. In the second 
and third algorithms kinetic energy is removed using the 
method described below. All three methods give equiv- 
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alent results for the total area. Unless otherwise noted, 
the results presented below were obtained with the third 
algorithm. We confirmed that the mean normal stress 
was independent of height to ensure that stress had equi- 
librated throughout the system. 

Within the Lagrangian framework, the finite-element 
discretization of the field equations leads to a discrete 
system of equations: 



final corrector configuration is therefore: 



Mx n+1+ F»'(x,x)=FJ 



(4) 



where x is the array of nodal coordinates, M the mass 
matrix, F^K the external force array, and FJJ+i the in- 
ternal force array arising from the current state of stress. 
The second-order accurate central difference scheme is 
used to discretize Eq. (0J in time [2^, [2(| . A small time 
step was used in order to be below the stability limit [25j . 

The above equations conserve energy and will not con- 
verge to the static equilibrium configuration. Optimum 
convergence is achieved by removing a fraction of the ki- 
netic energy of each node at regular intervals. The char- 
acteristic time for stress equilibration across the elastic 
cube is given by the time tj, for sound propagation across 
its height L. Equilibrium is reached in a few tl by scal- 
ing all velocities by a factor of 0.9 at intervals of tl/10. 
Other procedures gave equivalent results, but with longer 
run times. 

The internal forces F mt are calculated using a linear 
elastic isotropic constitutive law. All our results are ex- 
pressed in dimensionless form by normalizing pressures 
by the effective modulus E' . The Poisson ratio v was 
varied from to 0.45. Periodic boundary conditions are 
imposed at the contact surfaces to eliminate boundary 
effects. 

A contact algorithm is used only to enforce the impen- 
etrability constraint on the two surfaces. Adhesive and 
frictional forces are not considered in the current work. 
We adopt a conventional master/slave approach with 
a predictor/corrector split within the Newmark time- 
stepping algorithm 27] . The rough surface of the elas- 
tic cube is identified as slave while the rigid surface is 
master. The predictor part of the Newmark algorithm 
neglects the contact constraints and, therefore, consists 
of an unconstrained step, with the result: 



pred 
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(5) 
(6) 



This predictor solution needs to be corrected in order 
to comply with the impenetrability constraints. The 
net result of imposing these constraints is a set of self- 
equilibrated contact forces that modify the predictor po- 
sitions and velocities. Since the contact surfaces are pre- 
sumed smooth, normals are well-defined and the surfaces 
can be unambiguously classified as master and slave. The 
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(7) 
(8) 



Here denotes the nodal mass, the superscript Q S 
designates nodes that belong to the slave surface, and 
the vectors N and F are the normal and frictional forces, 
respectively. Friction will not be taken into account in 
the remainder of the paper, but is discussed in Ref. [27| . 

The formulation of an appropriate system of forces is 
obtained by considering a configuration in which a mas- 
ter surface triangle (facet of a tetrahedral finite element) 
is penetrated by several slave nodes. For each of the 
penetrating slave nodes, let 8 be the normal depth of 
penetration to be corrected by the contact forces. The 
contact constraints determine a local problem with the 
normal slave force as an unknown, which is obtained as 
a direct function of the penetration 6 and the master 
normal, n: 
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D. Calculating the Contact Area 



(9) 



The contact algorithm just described identifies all slave 
nodes on the top surface that attempt to penetrate the 
flat bottom surface. We obtain the total contact area by 
multiplying the number of penetrating nodes n p by the 
area of the square associated with each node. In most 
cases we report the fractional area A/Ao = n p /L 2 where 
L 2 is the total number of surface nodes. 

In principle, the area associated with each node varies 
with normal load if the Poisson ratio is non-zero. How- 
ever, tests using a Voronoi tessellation to determine the 
area for each node gave equivalent results for the rela- 
tively smooth surfaces considered here. A Voronoi ap- 
proach becomes important for rougher surfaces or irreg- 
ular grids. A more fundamental ambiguity in the contact 
area comes from the fact that contacts (Fig. |3J) contain 
many disjointed regions most of which contain only a few 
nodes. The consequences of this are discussed in Sec. 

There are also complications in defining the contact 
area in experimental systems. Dieterich and Kilgore's 
optical method identifies any region where the surfaces 
are closer than some fraction of the wavelength as in con- 
tact. Due to the fractal nature of contacts (Fig. |3J) this 
may overestimate the true area. At the atomic scale, con- 
tact is difficult to define. If one associates contact with a 
finite interaction, then A/A would always be unity since 
the van der Waals interactions between surfaces extend 
to arbitrary distances. A more practical definition is to 
associate contact with the separation at which the net in- 
teraction becomes strongly repulsive due to the overlap 
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FIG. 3: Regions of contact (dark) for a surface with L = 256, 
H = 1/2, v = 0, and A = 0.082. The fraction of the area in 
contact A/Ao = 0.1. 



0.2 p r-, r-^p p-i 
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0.002 0.004 0.006 0.008 0.01 



W/E'Aq 

FIG. 4: Fractional contact area A/Ao (solid line) as a func- 
tion of the normalized load W/E'A for L = 256, A = 0.082, 
v = and H = 1/2. The dashed line is a fit to the linear 
behavior at small areas. 



of electrons on opposing surfaces. This leads to a range 
of separations where surfaces are in contact. As in the 
optical measurements, the contact area can be greatly en- 
hanced relative to the penetration definition used here. 
Direct comparisons to atomistic models will be presented 
in another paper p8[ . 



III. RESULTS 

A. Contact area vs external load 

As noted in the introduction, analytic theories of con- 
tact between self-affine surfaces predict that the real area 
of contact A should b e propo rtional to the applied load 
W at small loads 0, El [l! . To make the load dimen- 
sionless we normalize it by Aq times the effective modulus 
E' , Figure0]shows a plot of the fraction of the projected 
area that is in contact A/Aq as a function of the normal- 
ized load for a system with L = 256, A = 0.082, v = 
and H — 1/2. As expected, the area is proportional to 
the load at small loads. A growing deviation from this 
proportionality is evident as A/Aq increases above 4 or 5 
%. 

To emphasize deviations from linearity, the dimension- 
less ratio of true contact area to load AE'/W is plotted 
against log 10 A/Aq in Fig. [5] Results for L between 
64 and 512 fall onto almost identical curves. The small 
variation between curves is comparable to that between 
different random surfaces of the same size. In each case 
AE'/W is nearly constant when from 1 to 8% of the sur- 
face is in contact. This implies that the mean pressure in 
contacting regions (p) = W/A is also constant. The ratio 



AE'/W drops as A/Aq increases further because A/Aq 
is bounded by unity while the normalized load keeps in- 
creasing. The ratio has roughly halved by A/Aq — 70%. 

Most of the analytic theories mentioned above explic- 
itly assume that there is a statistically significant number 
of asperities in contact, and that only the tops of asper- 
ities are in contact. The latter assumption breaks down 
as A/Aq approaches unity, contributing to the decrease 
at large A/Aq in Fig. [S] The first assumption must break 
down for our systems when the total number of nodes in 
contact, L 2 A/Aq, is small. This explains the rise in the 
L = 64 data for A/Aq < 2% in Fig. EJ This rise is de- 
pendent on the specific random surface generated, and 
is particularly dramatic for the case shown. Examina- 
tion of this and other data indicates that proportionality 
between load and area is observed when there are more 
than 100 contact ing nodes (L 2 A/A > 100 ). 

Batrouni et aL[l0j considered the same range of system 
sizes for v — 0. They fitted all data from A/Aq < 0.2 to a 
power law and found W cx A 1 with 7 > 1. Their numer- 
ical results clearly rule out an earlier prediction ^3 that 
7 = (1 + H)/2, but we do not believe that their results 
are inconsistent with an initially linear relation between 
load and area (7 = 1). Their fit included regions where 
the number of contacting nodes is below the minimum 
threshold just described. In addition, their value of 7 de- 
creased steadily towards unity with increasing L, varying 
from 1.18 at L = 32 to 1.08 at L = 256. We believe that 
this small difference from unity is within the systematic 
errors associated with the limited scaling range. Note 
that 7 = 1.1 would imply a 25% decrease in Fig. 0from 
A/Aq = 10~ 2 to 10 _1 which is much larger than the ob- 
served change (~ 5%). The inset of Fig. shows that 
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FIG. 5: The dimensionless ratio of the area to load AE' /W 
vs \og 10 (A/Ao) for the indicated system sizes. In all cases 
A = 0.082, v = and H = 1/2. The inset shows a linear plot 
of the same data. 



FIG. 6: The product k (Eq. II 01 as a function of roughness 
A for H = 1/2 and v — 0, and the constant values predicted 
by Bush et al. (solid line) [Tfll | and Persson (dashed line) @. 
The dotted line is a guide to the eye. 



data up to A/Ao = 0.7 can be described by a very sim- 
ple linear dependence of A/W on A. We conclude that 
7 = 1, but that leading quadratic corrections to scaling 
give a larger apparent exponent when a large range of 
data is fit. 

The lack of system-size dependence in Fig. may ap- 
pear surprising in the context of someprevious results 
for self-affine surfaces. These studies 0, QjJ considered 
a fixed roughness at large scales and examined changes 
in contact area with increasing resolution. They found 
that the contact area decreased as the number of nodes 
increased because the local slope of the surface became 
rougher at higher resolution. Our Fig. [SJcompares results 
for the same small scale roughness and shows that there is 
a well-defined thermodynamic limit as one increases the 
total system size. This result is not obvious, since the 
rms roughness at the scale of the entire contact rises as 
L H A. Apparently this increase in large scale roughness 
is irrelevant because it rises sufficiently slowly with L. As 
noted in Sec. Ill Al the large scale roughness is sensitive 
to the first few random numbers chosen in generating 
the self-affine surface. If surfaces with the same large 
scale roughness are compared, substantial differences are 
found because the small scale roughness varies. These 
fluctuations are absent when results from the same small 
scale roughness are compared. 

The existence of a well-defined thermodynamic limit 
allows us to consider results for a single system size in 
subsequent sections, and extrapolate the results to other 
cases. In the following sections we will focus on the con- 
stant region of Fig. [SJand examine variation of this value 
with the statistical properties of the interface. Unless 
noted, all results are for L = 256 and uncertainties due 



to statistical fluctuations are less than 5%. 



B. Comparison to analytic theories 

Bush et al. found that the ratio plotted in Fig. 
[3] should increase inversely with the root mean squared 
slope of the surface yj < \Vh\ 2 > . More specifically, they 
predicted that the quantity 

k = V< \Vh\ 2 >AE'/W (10) 

should have the constant value of \/2tt. Persson ar- 
rived at a rather different looking expression in terms 
of the height-height correlation function (Eq. 2] and |2J) 
0. Bowever, it can be reduced to a prediction that 
k = y^8/n using the fact that q 2 C(q) is the Fourier 
transform of |V/i| 2 . Note that both predictions have a 
well-defined thermodynamic limit that is independent of 
large scale roughness, just as observed in our simulations. 

For our surfaces, A is the rms change in height between 
adjacent nodes in each of the two spatial directions. Thus 
a/ < | V/i| 2 > = V2A. Numerical results for k are plotted 
against A in Fig. [f)|for v = and H = 1/2. The value of 
k only changes about 10%, while the roughness changes 
by almost an order of magnitude. If A is increased to 
larger values than considered here, the local slope of the 
surface exceeds unity in some regions. This regime was 
not studied because it requires a different meshing algo- 
rithm and most treatments of contacts assume that the 
local slope remains less than unity. 

We also evaluated k over the range of roughness expo- 
nents typically observed on real surfaces, 0.3 < H < 0.7 
0,0, and at the higher values of 0.8 and 0.9. Figure 
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FIG. 7: The product k (Eq. EU as a function of H for 
A = 0.082, and the constant values predicted by Bush et al. 
(solid line) [lR^ and Persson (dashed line) Q- The dotted line 
is a guide to the eye. 
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little effect on k. Results for each value of H are nor- 
malized by the value k(0) obtained at v = 0. In every 
case there is a nearly linear rise in k at small that 
appears to saturate as v approaches the limiting value 
of 0.5. The total change of around 10% is comparable 
to the change found with A. The increase in k with v 
appears to be related to increased interactions between 
nearby asperities. The lateral expansion in response to 
a normal stress increases with v. This reduces the local 
curvature, making it easier for adjacent regions to come 
into contact. Detailed analysis of neighboring asperities 
shows that a smaller peak between two contacts may be 
brought up in to contact at high Poisson ratios. 

All of the values of n in Figs. HO and lie between the 
analytic predictions of Bush et al. and Persson. Our re- 
sults suggest that using a value of k = 2.2 should predict 
the ratio of area to load within about 10% over a wide 
range of surface geometries at v = 0. Fig. [5] indicates 
that the value of k should be increased linearly to about 
2.5 as v increases to the limiting value of 0.5. 

The agreement with these analytic predictions is quite 
good considering the ambiguities in discretization of the 
surface. Both analytic models assume that the surface 
has continuous derivatives below the small length scale 
cutoff of the roughness. Bush et al. consider con- 
tact between elliptical asperities and Persson removes 
all Fourier content above some wavevector. While we 
use quadratic shape functions, the contact algorithm 
only considers nodal heights and assumes that contact 
of a node implies contact over the entire corresponding 
square. One might expect that this assumption would 
lead to larger areas of contact, and our results do lie 
above those of Persson. Discretization would not be im- 
portant if the spacing between nodes were much smaller 
than the typical size of asperity contacts. However, as 
we now show, the majority of the contact area consists 
of clusters containing only a few nodes. The number of 
large clusters grows as H — ► 1, which may explain why 
our numerical results approach Persson's prediction in 
this limit. 



FIG. 8: The ratio of k to its value at v = as a function of 
v. Results for H between 0.3 and 0.7 show nearly the same 
linear rise with v. Lines are linear fits to the data. 



[7| shows results for a fixed value of A = 0.082. Even 
though the large scale roughness, L H A, varies by more 
than an order of magnitude, n changes by less than 30%. 
There is a nearly linear decrease with increasing H that 
may be influenced by uncertainties in determining the 
true contact area. As shown in Sec. IIII CI increasing H 
also increases the population of large clusters and may 
reduce uncertainties associated with assuming that the 
entire square region around contacting nodes is in con- 
tact. 

Fig. [8] shows that the Poisson ratio also has relatively 



C. Distribution of connected contact regions 

Most continuum theories approximate the real contact 
area by summing over many disconnected asperity con- 
tacts, each of which has a circular or elliptical [l^ 
shape. The connected regions in our calculated contacts 
(e.g. Fig. I3J) are considerably more complicated. We 
consider two nodes to be connected if they are nearest 
neighbors on the square lattice of interfacial nodes [2^| . 
All clusters of connected contacting nodes are then iden- 
tified for each load. The area of each cluster a c is just 
the number of connected nodes, since each represents a 
square region of unit area. 

Figure El shows the probability of finding a cluster of 
a given area P(a c ) as a function of a c for H = 0.5, 
A = 0.082 and v = 0. Results for different system sizes 



8 



10° 

io- 1 

lO- 2 
IO" 3 

io- 4 

IO" 5 
10" 6 



1 1 1 1 1 II 1 1 


\ — i i i iii 




H=0.5 1 






f n 


i 


Cl. 
Iff- 


: 




- 

: 




|jK o -2 - 


_ 




r □ L=64 


\ ^ 


; a L=128 


°s : 

X 


r ■ L=256 




; ° L=512 
. i i 


i i ii.ii 



10 l 



10' 

a c 



10" 



10° 
10-11- 

io- 2 C- 

3 io- 3 t 

CO 

oT 10 -4g. 

10" 5 
10 -6 

10-7 



10° 



H=0.3 
H=0.5 J 
^- H-0.7 
-- H=0.9 4 




* b 



10 1 10 2 



10- 



FIG. 9: Probability P of a connected cluster of area a c as 
a function of a c for v = 0, H = 1/2, A = 0.082 and the 
indicated system sizes. All results follow a power law, P(a c ) ~ 
a,7 T , with r = 3.1 (dashed) line at large a c . The dotted line 
corresponds to r = 2. 



FIG. 10: Probability P of a connected cluster as a function 
of area a c for ^ = 0, A = 0.082, L = 512 and the indicated 
values of H. Dashed lines indicate the asymptotic power law 
behavior with r = 4.2, 3.1, and 2.3 for H = 0.3, 0.5 and 0.7, 
respectively. The solid line corresponds to r = 2. 



collapse onto a common curve. For a c > 8 the curve can 
be described by a power law P(a c ) ~ a~ T with t = 3.1 
(dashed line). This rapid fall off (t > 1) means that the 
integral of P is dominated by small clusters. Thus even 
though the maximum observed cluster size grows with L, 
the value of P at small a c is unaffected. All of the data 
shown in Fig. EJare for A/A between 5 and 10%, but we 
find that the distribution of clusters is nearly constant 
for A/Aq < 10%. This is the same range where A and 
load are nearly linearly related. The probability of large 
clusters rises markedly for A/Aq > 0.3, as clusters begin 
to merge and eventually percolate across the interface. 
The following data is all for A/A < 0.1. 

The model considered by Greenwood and Williamson 
also gives a load independent P(a c ) and this is a central 
reason for the linear relation between load and area in 
their model. As the load increases, each existing cluster 
grows larger and new small clusters are generated in a 
way that maintains a stationary distribution of cluster 
sizes. Only the total number of clusters changes, and 
it rises linearly with load. Our calculated P{a c ) and 
the total number of clusters both follow this behavior. 
However, the distribution and the shapes of the clus- 
ters are very different than assumed by Greenwood and 
Williamson. 

The variation of P(a c ) with H is shown in Fig. ITU1 
Results for very small clusters (a c < 8) are nearly inde- 
pendent of H, but the asymptotic power law behavior 
at large a c changes dramatically. The range of scaling 
behavior is too small for precise determination of r, but 
our data is consistent with r = 3.1 ± 0.2 for H = 0.5. 
When H < 0.5 there is an anticorrelation between the 
surface slopes in nearby regions [l^. This leads to more 



rapid up and down fluctuations that make large contacts 
unlikely and yield a larger r = 4.2 ± 0.4 for H = 0.3. 
When H > 0.5 there is a positive correlation between lo- 
cal slopes, yielding larger clusters. Fits give r = 2.3 ±0.2 
for H = 0.7 and r = 1.9 ± 0.1 for H = 0.9. For H < 0.9 
we find r > 2, which implies that the mean cluster size 
(a c ) is dominated by small clusters. Directly calculated 
sizes are indeed independent of both L and A/Ao. We 
find (a c ) = 1.8, 2.5, and 4.0 for H = 0.3, 0.5 and 0.7, re- 
spectively. For H = 0.9, r appears to be slightly smaller 
than 2, and the mean cluster size grows weakly with L. 

Some approximate treatments of contact begin by as- 
suming that the two surfaces do not deform and then 
determine the regions where the two solids would inter- 
penetrate 0, 0j El • Scaling arguments 0] and simula- 
tions 0, |3C| show that this rigid overlap model gives 
a power law distribution of connected areas at large a c 
with r = 2 — H/2. This is qualitatively consistent with 
Dieterich and Kilgore's experiments where r varied from 
1 to 2 and tended to decrease with increasing H . How- 
ever, it is qualitatively different from our results where 
t is always greater than 2. One consequence is that the 
mean cluster size from the rigid overlap model diverges 
with system size as L H , while it remains of order the 
discretization size in our calculation 29]. 

A possible explanation for the discrepancy between our 
results and experiment is that the latter identifies regions 
that are within some small fraction of the wavelength of 
light as being in contact. Figure ITU illustrates the dra- 
matic effect that this can have on P(a c ). The uppermost 
curve shows the cluster distribution obtained by apply- 
ing the rigid overlap model to our surfaces. The asymp- 
totic slope is consistent with the analytic prediction for 
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FIG. 11: Probability P of a connected cluster as a function of 
area a c for u = 0, H = 1/2, A = 0.082, L = 512 and different 
criteria for contact. The probability distribution for the rigid 
overlap model (open squares) falls off more slowly than a^ 2 
(dotted line). When only contacting nodes are included (h c — 
0), P falls off more rapidly (open circles). As the width h c 
of the region considered in contact increases, results from the 
full calculation approach the overlap results. 



H = 0.5: t = 2 — H/2 = 1.75. The lowermost curve is 
our result for the actual contact area. The intermediate 
curves were obtained by changing our definition of con- 
tact to include all nodes that are separated by less than 
some value h c . As h c increases, the number of nodes 
in the contact region rises and the probability of large 
clusters grows. When h c is comparable to or larger than 
A = 0.082, P(a c ) follows the rigid overlap prediction 
quite closely. It is likely that the optical experiments 
were in this limit. However, it is also likely that plas- 
tic deformation is important in these experiments. This 
effect will be explored in future work. 
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FIG. 12: Contact morphology for L = 256 under two different 
values of the displacement d following contact: d — 1.4, 2.54. 
Panels (a,d) show regions in true contact (h c = 0), while 
(b,e) show regions where the surface separation is less than 
h c — 0.1. Panels (c,f) show the contacts predicted by the 
rigid overlap model. The fractional contact area (A/Ao) is 
indicated for each case. 



D. Contact Morphology 

The contact morphologies produced by different mod- 
els are contrasted in Fig. ^| Results for two values of the 
interpenetration d are shown, where d is the downward 
displacement applied to the top of the elastic solid (Fig. 
|2J) after the surfaces first touch. The top panels show 
the results of the full calculation, the middle panels are 
for h c = 0.1, and the bottom panels show the contacts 
obtained from the rigid overlap model. The first obvious 
difference between the results is that the rigid overlap 
grossly overestimates the fraction of area in contact. For 
small d the actual area is roughly 8 times smaller than 
given by the overlap model. 

As seen above, the rigid overlap model also gives many 
more large clusters. Moreover the shapes of the clusters 



are quite different. Analytic studies predict that the over- 
lap model should give clusters with a non-fractal interior, 
but with a fractal interface. More specifically, if R is the 
diameter of a cluster, then the area a c oc R 2 , but the 
perimeter length s c o c R Df where the fractal dimension 
D f = (3 - H)/2 [lHH. Thus as the clusters grow in 
size, the perimeter becomes smaller and smaller relative 
to the area: s c /a c oc a c ( 1+H ^ 4 . Figure H31 shows that 
our results for the overlap model are consistent with this 
scaling prediction. However, the results for the full calcu- 
lation are quite different. The value of s c /a c approaches 
a constant at large a c indicating that the perimeter and 
area are both fractals with the same fractal dimension. 
Plots of s c and a c as a function of R show Df is roughly 
1.6 for H = 0.5. Note that s c is actually larger than a c 
because it is defined as the number of missing nearest- 
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to A are considered in contact, the distribution of clus- 
ters approaches that for the rigid overlap model. Panels 
(b) and (e) show the contact morphology produced in 
this limit (h c = 0.1). Note that the clusters from panels 
(a) and (d) have been connected into larger clusters that 
still lie within those of panels (c) and (f). Growth is most 
pronounced in regions where overlap is greatest. These 
regions carry a greater share of the load and flatten more. 

Greenwood and Wu p have recently reconsidered 
when a local peak should be considered as an asperity. 
They conclude that one should think of each cluster in 
the rigid overlap model as a single asperity, and use the 
diameter and height of the overlap to determine the di- 
mensions of an effective ellipsoidal asperity. Our results 
indicate that the original view [5j that almost all points 
are asperities and the typical asperity diameter is compa- 
rable to the lattice size provides a more accurate descrip- 
tion of the contact. However, the revised approach of 
identifying overlapping regions with asperities may give 
a better description of subsurface stresses, because it cap- 
tures correlations in the location of load bearing regions. 
Since the maximum shear stress is usually below the sur- 
face, the overlap model : 6] may be useful in modeling 
wear. 



neighbors along the periphery. Thus it would be 4 for a 
cluster containing a single node. 

Results for other values of H look very similar to those 
for H = 0.5. In each case s c /a c saturates at large a c , in- 
dicating the area and perimeter have the same fractal 
dimension. The limiting value of s c /a c decreases from 
about 1.7 for H = 0.3 to 1.1 for H = 0.9. The large 
number of perimeter nodes leads to some ambiguity in 
the total area obtained from our calculation. If only a 
fraction of the square around each node were actually 
in contact, then the true contact area would be smaller, 
movingour results in Figs. Hj] and [7] closer to Persson's 
result [1, |tj . It is interesting that our results approach 
Persson's prediction as H — > 1, and the perimeter be- 
comes less important. 

Despite the above differences, the rigid overlap model 
does provide information about where real contacts may 
occur. The distance between surfaces is always larger 
than that given by the overlap model because v = 0. 
Thus all of the contacts in panels (a) and (d) are part 
of the overlapping regions shown in (c) and (f). Only a 
fraction of the overlapping regions is in contact, because 
a local peak can screen a neighboring valley from con- 
tact. As pointed out by Greenwood [5|, a large fraction 
of points are local maxima, so the average cluster size is 
comparable to the lattice resolution. On the other hand, 
a small local maximum can only screen a small local re- 
gion. Thus there tend to be many small contacts in the 
regions where overlap first occurs. These points lie in 
the middle of the large clusters in panels (c) and (f). A 
higher density of clusters, and clusters of larger size, are 
found in these regions of panels (a) and (d). As noted 
above, when nodes that are within a distance comparable 



E. Distribution of local pressures 

Plastic deformation at the interface will be influenced 
by how pressure is distributed within the contact area. 
We find that this distribution has a strikingly universal 
form. Figure El shows that the probability P{p) for a 
contacting node to have local pressure p is independent 
of system size. Since the contact area increases linearly 
with load, the mean local contact pressure (p) ~ W/A is 
independent of contact area, and the entire distribution 
also remains unchanged for A/Ao between about .01 and 
.1. 

Increasing the small scale roughness (A) leads to a 
proportional increase in (p) . Yet Figure El shows that 
results for all A and H collapse onto a universal func- 
tion of the dimensionless variable p/(p). The proba- 
bility decreases monotonically with increasing p, and 
for p/(p) > 3 follows an exponential decay (solid line), 
P{p) oc exp(— p/pi), with pi ps (p)/1.6. This exponential 
tail implies that some regions have stresses much higher 
than (p) and may undergo plastic deformation even when 
the mean stress is much less than the hardness. 

Similar universal curves have been found for the stress 
distribution in a variety of "jammed" systems 123, in- 
cluding granular media [2(J, |2l| , thermal glasses [3l[ and 
polymer crazes In each case the tail of the dis- 

tribution follows a simple exponential rather than the 
Gaussian that might be expected from equilibrium ar- 
guments. Several explanations for the exponential form 
have been proposed [2(], 0, I^H |32, but most do not 
apply to our zero temperature, deterministic and elas- 
tic system. However, it is possible that the power law 
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FIG. 14: Probability distribution of the local pressure at 
contacting nodes for different system sizes with A = 0.082, 
v = 0, H = 0.5, and A/A between 5 and 10 %. 



correlations in interface height may lead to a hierarchical 
distribution of load that is analogous to the q- model [2l] . 

The distribution of local pressures plays a central role 
in Persson's theory of contact between self-affine surfaces 
8, 9J. He defines a resolution ( corresponding to the 
number of points along an axis at which the height of the 
surface is known, and assumes a smooth interpolation be- 
tween these points. Increasing £ corresponds to resolving 
more of the surface roughness and increases A if the sur- 
face is self-affine. The pressure distribution P(p, () is a 
function of both resolution and pressure. Its derivatives 
satisfy 



(11) 



where po = W/Ao is the apparent mean pressure, primes 
denote a derivative and 



G(C) = (£7'/Po) 2 <|Vh| : 



(12) 



Persson obtained solutions for P(p, C) in the geometry 
considered here by starting from perfectly flat planes with 
P(p, 1) = S(jj — pq) and iterating to higher resolution. He 
also imposed the boundary condition that the probability 
goes to zero at zero pressure |9(- 

Given the results shown in Fig. 15, it is interesting 
to ask if Persson's equations have a universal solution in 
the limit of small loads. Fig. 15 shows the probability 
distribution within the contact, while Persson's P(p, C) 
includes noncontacting regions and its integral over pres- 
sure is the fractional contact area. Thus it should be 
related to the universal distribution P by: P(p, £) = 
P{p/ < p >)p / < p > 2 , where p / < p >= A/A . 
Then use of Eqs. (10) and (11) leads to an equation for 
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FIG. 15: Probability distributions for p/(p) at the indicated 
values of A and H all collapse onto a universal curve. Here 
v = and A/Ao is between 5 and 10%. The solid line is a 
fit to the exponential tail of the distribution, the dotted line 
shows Eq. (14), and the dashed line shows a Gaussian with 
the appropriate normalization and mean. 



P: 



2P(x) + xP' (x) + k 2 P" (x) /4 = 0. 



(13) 



There is a solution with unit norm and mean for Persson's 
value of k — \/8/ir: 



~ 7T 7T 2 \ 

P{x) = -xexp(--a; ). 



(14) 



As shown in Fig. 15, this solution (dotted line) is much 
more strongly peaked than our numerical results, decay- 
ing to zero linearly in the limit p — > and as a Gaussian 
at large p. Note that since the mean and norm are the 
same for all curves, the presence of extra weight at large 
p implies more weight at low p. Also shown in Fig. 15 
is a pure Gaussian with unit norm and mean (dashed 
line). This solution provides a better fit to the numer- 
ical data at low p. However, as for "jammed" systems 
HHHHllllllll, the tail of the numerical distribution 
is much closer to a pure exponential than a Gaussian de- 
cay. This may reflect correlations in the loads carried by 
different asperities that are not fully captured in analytic 
theories 



asperities tn 

ilium 



IV. SUMMARY AND CONCLUSIONS 

In this paper, we developed a numerical framework 
for analyzing frictionlcss, non-adhesive contacts between 
self-affine surfaces using the finite element method. This 
method has been applied to perfectly elastic contacts 
with a range of Poisson ratios, roughness amplitudes and 
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roughness exponents. In each case the real contact area 
A rises linearly with load W until the fraction of the 
total area in contact reaches 5 to 10% (Fig. EJ. This 
implies that the average local pressure in the contacts, 
(p) = W/A, remains constant. The dimcnsionless pres- 
sure (p) I E' is independent of the system size even though 
the large scale roughness grows as L H . 

As predicted by analytic studies 0^|, the dimension- 
less pressure scales linearly with the small scale rough- 
ness A. The constant factor n (Eq. I1U|) that relates the 
roughness to the local pressure is always between the pre- 
dictions of Bush et al. and Persson Q (Figs. EHHJ- A 
value of k = 2.2 reproduces the numerical results within 
about 10% for v = 0, and the best fit value rises lin- 
early to about 2.5 as v rises to the limiting value of 0.5. 
These results allow the mean pressure and fractional con- 
tact area to be predicted for any elastic self-afnne surface 
with known small scale roughness. 

The detailed morphology of the contact region, and 
distribution of the areas a c of connected regions were also 
studied. As in early theories of contact [3, 113 j the in- 
crease in area with load reflects a linear increase in the to- 
tal number of contacts with no change in the probability 
distribution of contact areas. As W increases, each exist- 
ing contact grows, and new contacts are formed at a rate 
that maintains a constant P(a c ). At large a c , the prob- 
ability distribution falls off as a power law: P(a c ) oc a~ T 
(Fig. 1101) . Since r > 2 for H < 0.9, the mean cluster area 
is independent of L and comparable to the resolution of 
the calculation. 

The above results for connected clusters are consis- 
tent with the conclusion that a large fraction of nodes 
on a self-afnne surface are local maxima that should be 
treated as asperities J5j . However, recent experimental 
0,0 and theoretical U papers have suggested a different 
view. They examine regions where undeformed surfaces 
would overlap and associate each with a contact. This 
model gives qualitatively different distributions of areas 
(Fig. Hip . The value of r is always less than 2, and the 
mean cluster area diverges as a power of system size. The 
geometry of the clusters is also very different. The rigid 
overlap model gives two-dimensional clusters with fractal 
perimeters, while the full calculation gives fractal cluster 
areas with the same fractal dimension as the perimeter 
(Fig. 113(1 . Including regions where the surfaces are sepa- 
rated by less than h c as part of the contact leads to dra- 
matic changes in the cluster distribution exponent r and 
total area. The results approach the overlap distribution 
when h c is comparable to the small scale roughness. Op- 
tical experiments will detect gaps that are much smaller 
than the wavelength as in contact and this may explain 



why small values of r are observed. 

Plastic deformation will occur when the local pressure 
in a contact exceeds the hardness of the material. The 
linear relation between mean pressure and small scale 
roughness (Eq. 110(1 can be used to estimate when this 
will happen. The largest experimental values oip/E' are 
of order 0.1 and are obtained in amorphous and nanocrys- 
tallinc materials. Thus Equation 1101 implies that defor- 
mation can only be elastic when y/ (\ Vh\ 2 } < 0.1k ~ 0.2. 
This condition is violated for many surfaces, and the 
much smaller hardness of macroscopic crystals will lead 
to even tighter constraints on the roughness. Our ap- 
proach is readily extended to include plastic deformation, 
which will be the subject of future work. 

Plastic deformation may occur well before the mean 
pressure reaches the hardness because some nodes have 
local pressures much larger than (p). Results for all pa- 
rameters collapse on to a universal probability distri- 
bution P(p/(p)) (Fig. 115(1 . Persson has presented ap- 
proximate analytic equations for the pressure distribu- 
tion. This analytic distribution drops as a Gaussian at 
large p, while the numerical results have an exponential 
tail that greatly increases the number of sites with large 
pressures. Similar exponential distributions are found 
in man y ja mmed systems such as sandpiles, glasses and 
crazes [2(j, EI H3, El H3 • A common feature of these 
systems is a highly nonuniform distribution of stress. It 
is possible that the presence of small bumps on bigger 
bumps on still bigger bumps in our systems leads to 
stress transmission like that in the g-model for sandpiles 
pol l2l| . This hierarchical structure may produce stress 
correlations that are not included in the analytic model 


There are many interesting avenues for future research. 
The approach outlined here is readily extensible to in- 
clude more complex surface morphologies, plasticity, in- 
terfacial friction and tangential loading of the solids. 
More challenging issues include adhesion and the role of 
atomic scale roughness. These issues will require hybrid 
algorithms that include atomic information about inter- 
facial interactions. 
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